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ABSTRACT 

Context. Since 1997, BL Lacertae has undergone a phase of high optical activity, with the occurrence of several prominent outbursts. 
Starting from 1999, the Whole Earth Blazar Telescope (WEBT) consortium has organized various multifrequency campaigns on this 
blazar, collecting tens of thousands of data points. One of the main issues in the analysis of this huge dataset has been the study of 
colour variability. 

Aims. The massive amount of optical and near-infrared data collected during the campaigns enables us to perform a deep analysis of 
multiband data, with the aim of understanding the flux variability mechanisms. 

Methods. We use a new approach for the analysis of these data, focusing on the source spectral evolution. 

Results. We show that the overall behaviour of the BL Lacertae light and colour curves can be explained in terms of changing viewing 
angle of a moving, discrete emitting region, which causes variable Doppler boosting of the corresponding radiation. 
Conclusions. A fractal helical structure is suggested to be at the origin of the different time scales of variability. 

Key words, galaxies: active - galaxies: BL Lacertae objects: general - galaxies: BL Lacertae objects: individual: BL Lacertae - 
galaxies: jets - galaxies: quasars: general 



1. Introduction 



BL Lacertae (z = 0.0688 ± 0.0002; iMiller & Hawlevlll977l) is 
the prototype of a class of active galactic nuclei (AGNs), the BL 
Lac objects, which are well known for their pronounced vari- 
abilit y at all wavelengt hs, from the radio to the y-ray band (see 
e.g. I Villata et al.l I2002] for an extended description of this ob- 
ject). It has been one of the favourite targets of the Whole Earth 
Blazar Telescope (WEBTfl which has organized several multi- 
frequency campaigns on this source. 

The campaign carried out in May 2000 - January 2001 was 
particularly success ful, with periods of exceptionally dense sam- 
pling that allowed Villa ta et all (120021) to follow the intranight 
flux variations in detail. The slope of the optical spectrum was 
found to be weakly sensitive to the long-term brightness trend, 
but to strictly follow the short-term flux variations, becoming 
bluer when brighter. The authors suggested that the nearly achro- 
matic modulation of the flux base level on long time scales is due 
to a variation of the relativistic Doppler beaming factor, and that 
this variation is likely due to a change of the viewing angle. 

A subsequent campaign was carried out in 2001-2003. On 
that occasion, historical optical and radio light curves were re- 
constructed back to 1968, and both colour and time series anal- 
yses were performed. Colour an alysis on a longer time pe- 
riod c onfirmed the conclus ions by IVillata et al.l d2002l) . and al- 
lowed IVillata et al.l (l2004al) to quantify the degree of chroma- 
tism of both the short-term and long-term variability compo- 
nents. The optical spectral behavio ur was further analysed by 
Papadakis . Villata. & Raiteri (120071) . who interpreted the bluer- 
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when-brighter mild chromatism of the long-term variations in 
terms of Doppler factor variations due to changes in the view- 
ing angle of a curved and inhomogeneous jet. The study of 
the optical and radi o flux behaviour in the last forty years led 
Villa ta et al.l (|2009a) to propose a more complex scenario, where 
the emitting plasm a is flowing along a rotating helical path in a 
curved jet. Indeed, Rait eri etakl (12009) claimed that synchrotron 
plus self-Compton emission from a helical jet can explain the 
broad-band spectral energy distribution (SED) of BL Lacertae, 
from the radio to the y-ray band. 

We notice that in all the above papers dealing with colour 
analysis, the host galaxy contribution was subtracted from the 
total flux densities to analyse the colour behaviour of the AGN 
alone. 

A different approach was adopted by H agen-Thorn and coau- 
thors (see, e.g. Hag en-Thorn et all l2004t) . It is based on the as- 
sumption that the flux changes within some time interval are due 
to a single variable source. If the variability is caused only by its 
flux variations and the relative SED remains unchanged, then in 
the n-dimensional flux density space {F\ , F „} (n is the number 
of spectral bands used in multicolour observations) the observa- 
tional points must lie on straight lines. The slopes of these lines 
are the flux density ratios for the various pairs of bands. With 
some limitations, the opposite is also true: a linear relation be- 
tween observed flux densities at two different wavelengths dur- 
ing some period of flux variability implies that the slope (flux ra- 
tio) does not change. Such a relation for several bands would in- 
dicate that the relative SED of the variable source remains steady 
and can be derived from the slopes of the lines. 

That kind of analys is, using BVRIJHK phot ometry obtained 
during 1999-2001, led lHagen-Thorn et al.l d2004l) to the conclu- 
sion that the spectrum of the variable source in BL Lacertae can 
be described as F v oc v~ a , with different values of a for the op- 
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point out that this data set does not allow to discriminate whether 
a single variable source acts over the whole wavelength range. 

Th e lack of spectral variability found by Hag en-Thorn et all 
(2004) appears in contradiction with the bluer-when-brighter 
trend obtained in the previously mentioned works based on the 
WEBT data. We notice that a spectral flattening w ith brightness 
increase had already been found by Web b et all d!998l) when 
analysing BVRI observations of BL Lacertae around its 1997 
outburst. These authors suggested three possible explanations 
for the source spectral variability, involving changes in the elec- 
tron energy distribution, its turnover frequency, or the presence 
of multiple synchrotron contributions, but they were not able to 
favour one of them because of the lack of information at other 
wavelengths. A bluer- when-brighter behaviour was also detected 
by other authors in datasets more or l ess extended in time (e.g. 
IZhang et alll2004l:IStalin etafl 120061: lGulFaIll2006h . 

In this paper we investigate the matter with a larger dataset, 
to shed light on the source spectral variability and its possible 
origin. Our dataset is presented in Sect. 2, where light curves 
and flux-flux plots are shown and spectral variability discussed. 
In Sect. 3 we address the question of the origin of the variabil- 
ity, investigating the possibility that it is due to changes of the 
Doppler factor, which in turn are caused by orientation effects. 
Sect. 4 contains a summary of the results. 

2. Observational data and analysis 

One of the main goals of all the multiwavelength campaigns car- 
ried on during the last decades was the investigation of colour 
variability on different portions of the light curve, as well as 
the determination of correlations and temporal delays between 
variations at different frequencies. Apparently, the former issue 
can be easily addressed by acquiring good-quality photomet- 
ric data at different brightness levels, and then calculating the 
corresponding colour indices. However, the interpretation of the 
colour behaviour of blazars, particularly of BL Lacertae, still re- 
mains a topic of discussion by several teams of observers. It may 
be, at least partly, explained by the complexity of the source it- 
self, since one needs to subtract the radiation of the underlying 
elliptical galaxy, which is bright enough to contaminate the ac- 
tive nucleus photometry, especially at low flux levels. 

The optical and near-infrared light curves discussed in the 
following are shown in Figs.Q]and[2] They include all the WEBT 
data, complemented by Crimean and St. Petersburg optical data, 
and Campo Imperatore near-infrared data, acquired in 2005- 
2006 and 2008, i.e. in the periods not covered by WEBT cam- 
paigns. 

The WEBT data had already been carefully assembled and 
"cleaned". The colour analysis needs even more attention. We 
included for further processing only pairs of data that (1) were 
obtained with the same instrumentation within 15 minutes from 
each other; (2) did not show systematic differences from the 
most reliable data sets; (3) did not present scatter caused by low 
intrinsic accuracy. After this screening, we were left with the fol- 
lowing numbers of data pairs for each colour: 63 (U - R), 460 
(B - R), 1006 (V - R), 1004 (/ - R), 300 (J - R), 490 (H - J), 
and 473 (K - J). 

We transformed magnitudes into de-reddened flux den- 
sities using the Galactic absorption Ab — 1?42 from 
ISchlegel. Finkbeiner. & DavTsI d 19981). fh e lCardelli etail {l989) 



extinction law, and the lBessell et al.lJ 19981) flux -magnitude cali- 
brations. 

Figures|3]and|4]display flux-flux dependences for optical and 
NIR bands in the best-sampled time interval 2000-2008. They 
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Fig. 1. Optical light curves of BL Lacertae in 1994-2008. 



show that in no case the relation between flux densities can be 
properly fitted by a linear dependence, but rather by second- 
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order polynomials a + b ■ F R + c 
The concavity of the regression is always directed toward the 
higher-frequency band, thus clearly displaying a bluer-when- 
brighter type of variable sourc e behaviour. This res u lt is in agree- 
ment with thos e obtained by IVillata et al.l d2002l) . IVillata et al.l 
(l2004al) . and iPapadakis. Villata. & Raiteril (120071). bu t it is in 
contrast with what found bv lHagen-Thorn et al. (2004). Indeed, 
the latter authors used a limited dataset, spanning a smaller flux 
range, so they could not identify deviations from a linear depen- 
dence in the flux-flux plots. 

We assume that the deviations from regular variability seen 
in Figs. |3] and |4] are mostly due to noise. It cannot be excluded 
that some part of them are real, but their amplitude is small as 
compared to the overall variability pattern. 

Now we need to take into account the host-galaxy contri- 
butio n. Its flux densities are obtained by adopting R = 15.55 
from IScarpa et all (2000) and the mean colour indices for el- 
liptical galaxies from iMannucci et~aT] d2001l) . We get for the 
UBVRIJHK galaxy flux densities: 0.36, 1.30, 2.89, 4.23, 5.90, 
11.83, 13.97, and 10.62 mJy, respectively. The flux positions of 
the host galaxy are marked in Figs.|3]and|4]with diamonds. It is 
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Fig. 2. Near-infrared and 7?-band light curves of BL Lacertae in 
1999-2008. 



remarkable that these positions lie almost exactly on the extrapo- 
lations of the second-order regression curves for the correspond- 
ing pairs of bands. This suggests that whatever mechanism of 
variability acts, most likely it is one and the same for the whole 
optical-NIR region. 

All the measurements, according to WEBT prescriptions, 
were performed with a 8 arc sec aperture radius. Using a 
de Vacouleurs' radial flux distribution and typical values of see- 
ing, we estimate that * 50% of the galaxy flux density must be 
subtracted from the total photometric measure in order to get the 
flux density of the active nucleus. 

Fixing flux density values in R, the polynomial regressions 
of Figs.[3]and|4]allow us to get the corresponding flux densities 
in all the other bands. After subtracting the host galaxy contri- 
bution in each band, we finally obtain simultaneous spectra, free 
from the host galaxy and chance fluctuations. Figure [3] shows 
spectra of the variable source in BL Lacertae for four selected 
flux densities in R from 10 to 60 mJy. 



3. The origin of the variability 

In this section we investigate the possibility that the observed 
flux variations come from changes of the Doppler factor 5, and 
that these are in turn due to viewing angle variations of the emit- 



ting reg ion, as suggested by IVillata et al.l d2002l) . IVillata et"ai] 
d2004al) . and lViUata etaT1(l2009al) . 



We recall that 



6 = [T(l -jS-cosfl)]- 
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Fig. 3. De-reddened flux-flux dependences in the optical-NIR 
region. The red lines represent second-order polynomial regres- 
sions. Diamonds correspond to the position of the underlying 
galaxy. 
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Fig. 4. Same as Fig.[3]for the NIR region. 



/3 being the bulk velocity of the emitting plasma in units of the 
speed of light, F = (1 —0 2 )' 1 ' 2 the corresponding Lorentz factor, 
and 9 the angle between the velocity vector and the line of sight. 

In the observer's rest frame, the flux density is a func- 
tion of the Doppler factor 5 dRvbicki & Lightmaiil 119791 : 
lUrry & Padovanll 19951) 



FM=S p -F'M) 



(2) 



(1) 



where primed quantities refer to the source rest frame, and p = 
3 in case the radiation comes from a moving source (such as 
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Fig. 5. Spectra of BL Lacertae obtained by using the polyno- 
mial regressions shown in Figs.[3]and[4]and subtracting the host 
galaxy contribution. Numbers from 1 to 4 correspond to /?-band 
flux densities of 10, 20, 40, and 60 mJy, respectively. The pink 
lines tracing each spectrum represent cubic spline interpolations. 
Blue dashed arrows show the direction and amount of the spec- 
trum displacement when 6 changes by a factor 1.58 and p = 3 in 
Eq. 1 ; grey arrows correspond to p — 2 and a 5 variation of a fac- 
tor 1 .88. Blue dashed and black dot-dashed curves are the results 
of the shift of spectrum 4 toward spectrum 1 in both cases. 

a discrete, essentially point- like, emitting zone), while p — 2 
for a smooth, continuous jet (Begelma n et al.l[T9 84: Cawt hornel 
119911) . 

Doppler shift also affects frequencies: 

v = 6-V. (3) 

If we assume that the flux variations are due to changes in 
6, only in the case when the intrinsic spectrum is a power law 
F' V ,(V) oc (yy a and a is constant in the wavelength region under 
consideration, the flux-flux plots presented in the previous sec- 
tion would show linear dependences. Indeed, only in this case 
Doppler boosting effects on both F v and v lead to an unchanged 
observed spectral slope: 

F v (v) cc 5 p+a ■ v- a . (4) 

An example of that kind of behaviour was displayed by 3C 279 
in 2006-2007 dLarionov et al.L 120081) . However, in Sect. 2 we 
saw that the flux-flux relations are not linear, and indeed the 
spectra shown in Fig.[5]are far from a power law. 

As a consequence of Eqs. [2] and [3] in logarithmic scale (as 
in Fig. [5]l changes of A log 6 in the Doppler factor lead to a shift 
of the spectrum by A log v = A log 5 (in the frequency direction) 
and A log F v — p ■ A log 6 (in the flux direction). 

The blue dashed spectrum in Fig. [5] has been obtained from 
spectrum 4 by decreasing 6 by a factor 1.58, and setting p = 3. 
It fairly reproduces the "real" spectrum 1 . We cannot obtain the 
same agreement using p — 2 (black dot-dashed curve). This sug- 
gests that the general variability pattern of BL Lacertae during 
2000-2008 in the optical and near-infrared part of the spectrum 
can be interpreted in terms of Doppler factor variations affect- 
ing the radiation coming from a moving, discrete, and steady- 
emitting zone(s), possibly travelling in a jet. 

Under the assumption that 5 changes are due to viewing an- 
gle {8) variations, we can try to estimate what 8 variations are 




e, degrees 



Fig. 6. Histogram of viewing angle 8 distribution in 2000-2008, 
fitted with (a) two Gaussians; (b) ring-like Gaussian distribution 
plus Gaussian (red dashed curves indicate the contributions of 
the different components); (c) ring-like distribution with asym- 
metric density (see text). 



needed to account for the observed range of flux variability in 
the R band. From the definition of 6 in Eq. [1] 

ST- 1 

8 — arc cos . (5) 

5Vr 2 -i 

With p = 3, Eq. |4]becomes F v (v) = F 5 3+a v~ a . The mean 
value of a close to the R band is a a 1.3 (see Fig. Q. The 
constant F can be determined from F a = F max v" /6^£, where 
^max is the maximum observed flux density in R band and (5 max 
the corresponding maximum Doppler factor, which comes from 
the minimum viewing a ngle 8 m \„. We fix, as a te ntative value, 
8 mia = 2°. According to Marscher et alj d2008l) and lJorstad et al.1 
(2005), r = 7.0 ± 1.8, then 6» min = 2° yields 6 max = 13.15. 

Using these values in the above equations, we can get the 
pertaining values of 8 for each observed 7?-band flux density. 
The histogram of the distribution of 8 for all the time interval 
(2000-2008) considered is shown in Fig. [6] In order to avoid 
false signals caused by uneven distribution of the data, we made 
nightly binning beforehand. 

Figure [6] shows that all the values of 8 are grouped within 
the range from 2° to 8° (notice that the value of 2° was fixed 
initially; setting another value would just shift all the histogram 
leftwards or rightwards). The resulting distribution shows sub- 
stantial asymmetry and bimodality. Of course, the interpretation 
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of these features is not unambiguous. Red dashed lines in the 
two upper panels represent fitting with (a) two Gaussians and 
(b) ring-like structure plus Gaussian; blue vertical lines corre- 
spond to the centres of these components and the solid black 
lines indicate their sum. The bottom panel (c) represents the 
single-component case of a ring-like distribution with increasing 
density toward larger angles. These three model distributions are 
shown in Fig. [7] 

In order to discriminate among them, we recall that we 
are working under the assumption that the flux variations of 
BL Lacertae are mainly (if not solely) due to changes of the 
Doppler boosting. In this case it is necessary to take into ac- 
count also time contraction, that leads to shortened time scales 
in the observer's rest frame: t = 6~ l t'. The degree of contrac- 
tion is inversely proportional to the Doppler factor, which de- 
pends on the viewing angle of the velocity vector of the emitting 
source. It means that the smaller is during some interval in the 
light history of BL Lac, the faster this interval will pass in the 
observer's rest frame. Quantitatively, in our case it implies that 
when 9 changes from 8° to 2°, the time scales shorten by a factor 
~2. 

This argument favours the model of Fig. [7b: the circular 
shape ma y reflect dominant sp iral movement of emitting features 
(see, e.g. lMarscher et alll2008l) . with time contraction leading to 
shorter life-time of high-flux-level events. This can also be seen 
in Fig. [8] where we plot the time evolution of 6 corresponding 
to the 1-day binned light curve in 1994-2008: smaller angles 
appear to last for a shorter time. 

An alternative explanation o f both Figs. [ 7 b and [8] can be 
found in the helical-je t model by Villata et al. (2009 al see also 
Villata & Raite rilll999h iRaiteri et al.ll 19991: lOstorero et al.ll2004t 
Raite ri et al. l2009Tr In this scenario, the changing-angle emit- 
ting zone is the whole optical-NIR emitting region in the heli- 
cal jet, whose orientation changes because of the helix rotation. 
In this case Fig. [7b would represent all the viewing angles ac- 
quired by the emitting zone during this rotation in 2000-2008. 
The angle time evolution in Fig. [8] is rather noisy and does not 
allow to recognise any clear periodicity, even if something like 
a yearly modulation can be envisaged (see also Fig. 1). The ori- 
gin of the noise around the main modulation may reflect a more 
complex structure of the helical path, such as a thinner helical 
structure of which the main helix would be the axifl In this pic- 
ture, fast flares in the light curve (i.e. the noise in the angle curve 
and the thickness of the ring in Fig. [7b ) could be due to shocks 
moving along the fine helical path and hence rapidly changing 
their viewing angles. When these events occur during a mini- 



2 This thinner helical structure may simply consist of helically 
twisted magnetic field lines, whose toroid al component is neede d for 
the MHD equilibrium of the jet (see e.g. IVillata & Ferraril 1 19951 , and 
references therein). 



Fig. 8. Time history of the viewing angle 6 from 1994 to 2008. 
The red line shows a cubic spline interpolation through the 30- 
day binned data. 



mum viewing angle orientation of the whole emitting region in 
the main helix, we observe the optical-NIR outbursts. 

From Figs. 1 and 8 one can also see tha t before 1997 the opti - 
cal activity was m uch lower. According tolVillata et alj d2004b), 
iBachet al.1 d2006l) . and lVillata etai1d2009al) the 1994-1996 faint 
state shown in the figures repr esents just the end of a longer qui- 
escent period started in 198 1. 1 Villata et al.l d2009af) ascribed this 
low-moderate activity to a general bad alignment of the opti- 
cally emitting region with the line of sight. Namely, the axis of 
the main helix had a viewing angle larger than that acquired in 
1997. (Looking at the figures, it seems that this angle has been 
again increasing after 2004.) According to this picture, if we had 
sufficient data in 1981-1996, we could compose a figure similar 
to Fig. [7b, but with the ring centre (representing the helix axis) 
farther from the line of sight. Thus, Fig. [7b would be a neces- 
sary simplification that does not take the helix drift into account, 
and this may be the reason why the single ring does not fit well 
the "wings" in Fig. 6c, which would belong to other, closer and 
fa rther, rings. 

IVillataetail d2009a) also argued that the helix axis is curved, 
to explain why the optical low activity was actually accompanied 
by strong low-frequency ra dio activity, and vice versa (see also 
IVillata el all 120071 l2009bl for the case of 3C 454.3). If the he- 
lix axis is curved, it may even be helical itself. In conclusion, 
we could be in the presence of a "fractal" helical structure: the 
finest structure being responsible for fast flares (on day-weeks 
time scales), the intermediate one for months-lasting, roughly 
annual outbursts, and the biggest structure for the longer-term 
(ten-yearly or more) alternation of low and high activity. 

4. Conclusions 

The analysis of the photometric U BVRIJHK data collected dur- 
ing the 2000-2008 WEBT campaigns on BL Lacertae, comple- 
mented by Crimean, St. Petersburg, and Campo Imperatore data 
in 2005-2006 and 2008, led us to the following conclusions: 

1 . The distribution of the optical and NIR data on flux-flux plots 
indicates the presence of a single dominant source of vari- 
ability with a bluer-when-brighter behaviour. 

2. The dependence of the spectral shape on the flux level agrees 
with the suggestion that the variability in that spectral region 
is mostly caused by changes of the Doppler factor. 

3. The emission would likely come from a moving, discrete, 
and steady-emitting zone(s), possibly travelling in a jet. 
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4. Under the assumption that the Doppler factor variations are 
due to changes of the viewing angle, the histogram of view- 
ing angle distribution appears asymmetric and can be fit- 
ted with a ring-shaped pattern with uneven angular distri- 
bution. We infer that this is caused by time contraction in the 
observer's rest frame, which sharpens the high-flux (small 
viewing angle) events. 

5. The noise around the main viewing angle modulation could 
be the signature of a finer jet structure. The above men- 
tioned discrete emitting zones (e.g. shock waves) would 
travel along a short-pitch helical path, thus rapidly chang- 
ing their viewing angle to produce fast flares. The axis of 
this helix would be helical itself, and its changing viewing 
angle due to rotation would cause the roughly annual out- 
bursts. This main helix, in turn, would have a curved axis, 
maybe helical again, whose orientation changes would yield 
the observed long-term alternation of high and low optical 
activity. In summary, the jet might exhibit a fractal helical 
structure, where increasing geometric scales correspond to 
increasing variability time scales. 

The analysis presented in this paper is based on data cover- 
ing a limited energy range (optical-NIR). However, our results 
confirm the geometrical interpretation of BL Lacertae variabil- 
ity suggested in previous works based on WEBT data, where 
broader-band observations were available. We cannot rule out 
other explanations for the observed spectral variability, based on 
intrinsic rather than geometric mechanisms. However, we think 
that straight and motionless jets are not realistic (as the bent tra- 
jectories of jet components revealed by VLBI observations sug- 
gest), and thus changes in the viewing angle of the emitting jet 
regions are very likely to occur. We expect that in the near future 
broad-band multiwavelength campaigns, including y-ray obser- 
vations by the Fermi-GST and AGILE satellites, will help clarify 
the picture on blazar variability. Indeed, they will allow people 
to study simultaneously both the low-energy synchrotron and the 
high-energy inverse-Compton emission components, and thus to 
better test theoretical models. 
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